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Abstract 

We develop a detailed theoretical framework for various types of transcription factor gene oscillators. We further 
demonstrate that one can build genetic-oscillators which are tunable and robust against perturbations in the critical control 
parameters by coupling two or more independent Goodwin-Griffith oscillators through either -OR- or -AND- type logic. 
Most of the coupled oscillators constructed in the literature so far seem to be of -OR- type. When there are transient 
perturbations in one of the -OR- type coupled-oscillators, then the overall period of the system remains constant (period- 
buffering) whereas in case of -AND- type coupling the overall period of the system moves towards the perturbed oscillator. 
Though there is a period-buffering, the amplitudes of oscillators coupled through -OR- type logic are more sensitive to 
perturbations in the parameters associated with the promoter state dynamics than -AND- type. Further analysis shows that 
the period of -AND- type coupled dual-feedback oscillators can be tuned without conceding on the amplitudes. Using these 
results we derive the basic design principles governing the robust and tunable synthetic gene oscillators without 
compromising on their amplitudes. 
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Introduction 

Transcription factors (TFs) regulate the quantitative levels of 
several proteins inside a living cell [1^]. TF networks present 
across various organisms ranging from prokaryotes to higher 
eukaryotes and consist of fundamental building blocks such as 
autoregulatory loops, cascades and single input modules, feed- 
forward and feedback loops, dense overlapping regulons and 
oscillatory loops [5-7] . Feedback loops act as bistable switches and 
feedforward loops have been shown to act as efTicient filters for 
transient external signals [8], [10-12]. Positive self-regulatory 
loops seem to play important roles in the maintenance of cellular 
memory [3] and subsequent reprogramming of the cellular states 
whereas negative auto regulatory loops have been shown [11] to 
speed up the response times against an external stimulus [8-10], 
[12]. Oscillatory loops drive the developmental as well as mitotic 
cell-cycle dynamics [13] and circadian-rhythms [14], [15] 
associated with the intracellular concentration of various types of 
proteins, metabolites and other cell-signaling molecules. Under- 
standing of the detailed dynamics of oscillatory loops associated 
with the TF networks is a central topic in biophysics, synthetic and 
systems biology. 

The minimalist TF network model that can generate self- 
sustained oscillations is the well-known Goodwin-Griffith oscillator 
which has a single gene that codes for a TF protein that negatively 
auto-regulates its own transcription [16-18]. In this model the TF 
protein-product undergoes a one-step modification that yields the 
matured or active end-product and subsequendy n numbers of this 
end-product bind with the cis-regulatory modules (CRMs) of the 
associated promoter that in turn results in down-regulation. Here 



n is the Hill coefficient associated with the cooperative type 
binding. Detailed studies on this minimalist model showed [17] 
that the inequality condition n>8 is necessary to generate self- 
sustained oscillations in the levels of mRNA and protein. This 
result was obtained with the assumptions that the rate constants 
associated with the synthesis and decay of the protein and end- 
product are equal and the binding-unbinding of the end-product 
with the promoter is much faster than the rate of change in the 
synthesis and degradation of mRNA, protein and end-product. 
Further it was assumed that the decay of mRNA and protein 
product follows a first order type reaction. 

It was realized later that the inequality condition n>8 is unlikely 
[19] under in vivo conditions since the formation of such large 
multimeric protein complexes via pure three dimensional diffirsion 
(3D) limited collisions (Figure 1) is almost an improbable event 
and several other modifications over the Goodwin-Griffith model 
were proposed to reduce the required value of n. Most of these 
modifications were mainly associated with the insertion of (a) a 
temporal delay in the negative auto-regulatory loop either 
exphcidy as a time-delay in between the synthesis and binding of 
end-product at the promoter [19], [20] or implicitly via inserting 
additional reaction steps [19] in the formation of end-product that 
interacts with its own promoter and (b) a non-linear Michaelis- 
Menten type kinetics in the decay of mRNA and protein products 
despite of the first order type kinetics and (c) additional TF gene 
members in the negative feedback loop which again indirectly acts 
as temporal delay in the overall negative feedback. The delayed 
negative feedback may also be coherently or incoherently 
amplified [21-23] via the insertion of a positively regulated 
intermediate. Here the temporal delay is connected with the 
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overall time that is required for the transport of fully transcribed 
mRNAs from nucleus to cytoplasm, post-translational modifica- 
tions and subsequent transport of active TF proteins into the 
nucleus through 3D diffusion. When the decay of mRNA and 
protein product follows a Michaelis-Menten type kinetics then the 
Goodwin-Griffith (GG) oscillator seems to produce self-sustained 
oscillations [19] even for n= I. The genetic oscillator module with 
three TF genes connected in a cyclic negative feedback loop is 
named as repressilator [24]. Though these motifs were shown to 
be oscillatory through deterministic and stochastic simulations, 
significant fraction of cells containing the constructs of these motifs 
were not showing any oscillations under in vivo experimental 
conditions. It was argued that it could be partially due to the noisy 
nature of intracellular environment [18], [24]. Here one should 
note that most of the simulation studies were performed with 
constant parameter values which may not be true under in vivo 
conditions. In this context it is essential to investigate how the 



oscillatory dynamics of these motifs reacts to perturbations in the 
system parameter values. 

Most of the earlier studies on GG and other oscillator models 
assumed a quasi-equilibrium condition for the binding-unbinding 
dynamics of the negatively autoregulated TF proteins at their own 
promoters. This is mainly to reduce the four or higher dimensional 
Jacobian matrix associated with the non-linear system of 
differential rate equations into a three dimensional one to ease 
further analysis since there is an additional rate equation 
corresponding to the promoter state dynamics apart from the 
rate equations associated with mRNA, protein and end-product. 
However this assumption is valid [8] , [9] only when the timescales 
associated with the synthesis and degradation of mRNAs and TF 
proteins are much slower than the timescales associated with the 
binding-unbinding of regulatory TFs at the respective promoters. 
Recent studies [8] on feedforward loops suggested that the 
binding-unbinding dynamics of TF protein at the promoter can be 
ignored only when the cellular volume Vc ( = volume of nucleus in 
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Figure 1. Goodwin-Griffith genetic oscillator model. The transcription factor (TF) gene A is transcribed and translated into the TF protein 
product that in turn is converted to the active end-product. The end-product (or its oligomer as in case of lad repressor negative-feedback-only 
system that was constructed in reference [25]) is the key molecule that locates the respective c/s-regulatory elements associated with the promoter of 
TF gene A through a combination of one-dimensional (1 D) and three-dimensional (3D) routes as that of typical site-specific DNA-protein interactions. 
Here either monomers of the end-product directly assemble at the corresponding regulatory elements in a combinatorial manner (I) or the fully- 
formed complex of n^-mer binds with the respective regulatory sites (II). Assembly of combinatorial TF molecules results in the looping of DNA 
segment that is present in between the promoter and c/s-regulatory elements. ARPC is the assembled repressor-promoter complex that in turn 
results in down-regulation. Our analysis shows that out of these two competing pathways, the pathway I is the most probable one since it takes 
shortest time. The corresponding set of differential equations is given in Eqs (4-5). This system is well characterized by parameters of Group I, II and 
III. Group I consists of parameters (u'„, >'„,£„) whereas Group II consists of equilibrium parameters (x„,;(„) and Group III consists of ordinary type 
perturbation parameters ((j„, ;<:„,/„). Most of the earlier studies assumed zero values to Group II parameters apart from assuming zero for v„ that 
controls the promoter state dynamics. Blue colored spheres are the dimers of lac repressor. Here a.a denotes amino acids and n.a denotes 
nucleotides. 

doi:1 0.1 371 /journal.pone.01 04328.g001 
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case of eukaryotes) is comparable with that of the prokaryotes [8] 
such as E. coli (F^ ~10 m^) and the influence of the promoter 
state fluctuations on the overaU dynamics of feedforward/feedback 
loops seems to significantly increase as the nuclear volumes 
increases as in eukaryotic cells across yeast to human. Further, the 
Michaelis-Menten type degradation kinetics associated with 
mRNA and protein is a valid assumption only when the 
concentrations of these species are much higher than the 
concentration of the corresponding nucleases and proteases. 
Nevertheless in most of the in vivo conditions, the intracellular 
levels of mRNA and protein of a particular TF gene will be much 
lesser than the corresponding levels of the non-specific nucleases 
and proteases. When the latter is true then the enzyme mediated 
decay of mRNA and protein will eventually follow a first order 
type kinetics. In this article, using a combination of theoretical and 
simulation tools (a) we develop a generalized theoretical frame- 
work of various types of genetic oscillators by explicitly incorpo- 
rating the promoter state dynamics and other chemical reaction 
balances in detail. Using our detailed model (b) we identify and 
classify various critical control parameters and compute their 
physiological ranges which are required to generate self-sustained 
oscillations in the intracellular levels of mRNAs and transcription 
factor proteins and (c) explore various possibilities of coupling 
independent gene oscillators and fine-tuning the period of such 
coupled system. We further (d) demonstrate that by coupling two 
or more independent Goodwin-Griffith oscillators one can design 
oscillatory network architectures which are tunable and also robust 
against perturbations in system parameters. 

Results 

Theoretical framework of transcription factor gene 
oscillators 

The Goodwin-Griffith oscillator consists of a negatively self- 
regulated gene (we denote it as TF gene A) which codes for a 
transcription factor protein (Figure 2A). We denote the cellular 
concentrations (mol/lit, M) of its mRNA as protein as p^, the 
transformed c-nd-product as and the complex of promoter with 
the end-product as Xa- Here the total ceUular concentration of 
promoter is d^ji and the overaU promoter state occupancy by the 
end-product will be X„ =x,Jd,,i where X„ e (0, 1). Though there is 
only one copy of the promoter by definition, we use a continuous 
type probability variable X„ to describe promoter state occupancy 
mainly to account for its partiaUy occupied status [8], [9]. The 
transcription and translation rates are denoted as k^a (Ms ') and 
kpa (s ') respectively. The first order decay rate constants 
corresponding to mRNA and TF protein are y„,a (s ') and yp„ 
(s ') respectively. The first order on- and ofF-rates associated with 
the transformation of protein into the matured end-product are 
denoted as kaf (s ') and kar (s ') and the corresponding 
dimensionless dissociation constant is ka = Xar/kaf. The overall 
forward and reverse rate constants associated with the binding and 
unbinding of w^, numbers of end-product molecules with the 
respective CM-regulatory modules (CRMs) of the promoter of TF 
gene A are k„f (M^"".s^') and kar (s ') and the corresponding 
dissociation constant is defined as Karf=kar/kaf (M""). To 
simplify the analysis further we introduce the foUowing scaling 
transformations to project the time and concentration variables 
into the dimensionless space. 



'^ = ypat\Pa=Pa/Pas; Ma=ma/mas; = zjpas; Xa= xj da^: (1) 

In these equations z denotes the dimensionless time that is 
measured as the number of lifetimes of the protein product of TF 
gene A and P„, M„, Z,^ and X,^ are respectively the dimensionless 
concentrations of protein, mRNA, end-product and promoter 
complexes. We also should note that (Pa, M^, andX^) e (0, 1) by 
definition. Here the steady state values of mRNA and protein in 
the absence of negative self-regulation can be defined as foUows 
[8], [9]. 

Pas = kmakpa / Ymaypa'^ f^as=kma/yma (2) 

We further transform the parameter associated with the 
multimerization of end-product and subsequent binding events 
as foUows. 

Ma = Karf/P2 Karf = kar j kaf (3) 

Using the scaling transformations given by Eqs (1-3) one can 
write the deterministic rate equations corresponding to the 
temporal evolution of dimensionless concentration variables (X^,, 
Ma, Pa, and over dimensionless time variable x as follows. 

VadXaldx = Za''-{\-Xa)-HaXa 

WndMn/dx = (l—X„)— M„ 

(4) 

dPa/dx = Ma-Pu- Oa{Pa - W 

eadZa/dT = Pa - {Xa + Ka)Za- XaiZa"" {I - Xa) - Ha^a) 

The initial conditions are {Xa,Za,Ma,Pa) = 0 at T = 0. When 
[Va = 0, (Ja = 0 and Ka = 0), then this system reduces to the usual 
GG oscillator model for three concentration variables. Here we 
have defined the dimensionless ordinary perturbation parameter 
Ka = y,aj).af where y-„ (s~') is the first order decay rate constant 
associated with the protein end-product Z^. Since y^^ » Xaf will be 
true in most of the physiological conditions and Ka is an ordinary 
perturbation parameter one can assume Ka~0. When there is a 
dimerization of z„ (we denote the dimer Zg-i,, as j,^ as in case of 
Lac r(;press()r system that has been constru('t(;d and studied in 
reference [25] (negative-feedback-only model using lad gene) then 
the first and last equations of Eqs (4) wiU be modified as foUows. 

VadXaldx = Ya"" ( 1 " X«) - /i„X„ 

S„dZa/dx = Pa- (4 + Ka)Za - ffya [ZaZa - Xya Y a) (5) 
^yadYaldx = ZaZa " lya Ya " lya{ Ya"" ( 1 " X,) - /X,X,) 

Here various parameters associated with the dimerization of 
end-product of TF protein A and subsequent assembly of this end- 
product at the own promoter are defined as foUows. 
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Figure 2. Various types of transcription factor based genetic oscillators. In case of positive regulation the combinatorial transcription 
factors bound at c/s-regulatory modules enhance the initiation of transcription by strengthening the RNAPII-promoter interactions through their 
distal action (positive arrows) whereas in case of negative regulation, the RNAPII-promoter complex will be destabilized by the combinatorial TFs 
present at CRMs (negative arrows). A. Goodwin-Griffith oscillator. Bl. One-to-one dual feedback oscillator. Here the end-product of TF gene A binds 
at the promoter of TF gene B and down-regulates it whereas the end-product of TF gene B binds with the promoter of TF gene A and down-regulates 
it. B2. Two independent Goodwin-Griffith oscillators are coupled through -OR- type logic with NN-NN configuration. Here the promoter of TF gene A 
will have binding sites for the end-products of both TF gene A and B and so on for TF gene B. 83. GG oscillators are coupled through -AND- type logic 
with N-N configuration. 84. GG oscillators are coupled through -OR- type logic with NN-PP configuration. 85. GG oscillators are coupled through -OR- 
type logic with NP-NP configuration. 86. Possible robust synthetic gene oscillator. Here K is the booster TF gene that is coupled to N-N type dual 
feedback oscillator via -OR- gate. CI. Repressilator that is built with three TF genes by cyclic coupling. C2. Three independent Goodwin-Griffith 
modules are coupled through -OR-type logic. Here dashed lines show the fully interconnected network. C3. Three independent Goodwin-Griffith 
modules with -AND-type logic. Here dashed lines show the fully interconnected network. 
doi:1 0.1 371 /journal.pone.01 04328.g002 



kya = /^rya I Pas^fya] ^ya —Pas^fya/ \fa 

Here we have defined Ya=yalPas and /./,.« (M~' s~') is the 
forward rate constant associated witli tlie dimerization reaction 
and Irya (s ) is the corresponding reverse rate constant. One 
should note that for a fully functional Lac repressor system the HiU 
coefficient will be ra^ = 4 since an octamer of Lac repressor protein 
(which is a dimer) is involved in the overall looping of DNA that 
results in strong repression of lad. The system of Eqs (4) is 



completely characterized by the following set of dimensionless 
parameters. 

'^a = 7pff/ Pas^af) H'a = 7 pa/ yma\ £fl = Jpa/ ^af\ 
da = Kfl Jpu ; Xa =Pm ' ' da-_k„f / Xaf 

Here one should note that the parameters {fia^Xa'^a) '^'"C 
functions of n„ that can be simplified by assuming an in vivo 
protein level as Pas = 1 • To simplify the analysis further we can 
classify these dimensionless control parameters into Group I, II 
and III. Group I consists of (va,M'a,Ea) which are all singular type 
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perturbation parameters since they multiply the first order 
derivative terms. One should note that this set of parameters 
directly controls the dynamics of changes in the cellular 
concentrations of active promoter, mRNA and end-product 
respectively. Group II consists of (ffa, k„) those are ordinary 
type perturbation parameters. In Group II, Gg controls the 
coupled dynamics associated with the concentrations of TF protein 
A and its end-product whereas controls the coupled dynamics of 
changes in the concentrations of end-product and it's binding with 
the promoter sequence. The lifetime of end-product is controlled 
by Ka- One should note that almost all the earlier studies assumed 
that (ffo, Xa; Ka) - 0. Group III consists of the equilibrium and 
promoter affinity parameters /ij. In this controls the 
equilibrium associated with the formation of end-product and 
controls the equilibrium associated with the binding of Ua 
molecules of end-product with CRMs of the promoter of TF 
gene A. 

Biophysical modeling of promoter state dynamics 

The total time required to initiate transcription consists of at 
least two different components viz. the time required (proportional 
to 1/A„y) for the assembly of n„ numbers of TFs at the respective 
CRMs of promoter and the time required for subsequent looping 
of DNA and subsequent distal action of TFs on RNAPII-promoter 
complex. The time component associated with the looping and 
distal action along with the time required for elongation and 
termination steps are included in the definition of transcription 
rate (the total time required for transcribing a full length mRNA 
win be equal to \lk^!T\ie. kinetics of interaction of molecules 
of end-product with the sequentially located CRMs can occur in 
two different possible ways namely binding of the full-fledged 
complex of Mo molecules of end-product (pathway II) or sequential 
assembly of the monomers of end-product at the corresponding 
cis-regulatory DNA-binding sites similar to that of a combinatorial 
binding of TFs with CRMs (pathway I) as in eukar^'otic systems 
(Figure 1). Though the pathway I resembles a (Wa+1)* order 
reaction it is stiU an operable one since the length scale of the 
genomic DNA is much higher than the combinatorial binding TF 
proteins. Binding of Wa numbers of transcription factors in a 
sequential manner or Ma-mer of end-product leads to looping of 
the DNA segment that is present in between promoter and CRMs 
of TF gene A that results in the spatial or distal communication 
between the end-product present at CRMs and the already formed 
RNAPII-promoter complex which in turn activates (positive) or 
deactivates (negative) the initiation of transcription depending on 
the type of self-regulation [1], [3], [4], [26], [27]. In case of 
activation or positive regulation, the combinatorial transcription 
factors bound at CRMs enhance the initiation of transcription by 
strengthening the RNAPII-promoter interactions through their 
distal action (positive arrows in Figure 2) whereas in case of 
negative regulation, the RNAPII-promoter complex will be 
destabilized by the combinatorial TFs present at CRMs (negative 
arrows in Figure 2). Here the destabUization of RNAPII- 
promoter complex may be through the formation of stem and 
loop structures. In prokaryotes, these types of up and down 
regulations generally do not involve recruitment or combinatorial 
binding of several TFs and the regulator transcription factor 
directiy influences the RNAP-promoter interactions as in case of 
negatively self-regulated oscillatory motifs constructed with a lac- 
repressor gene. Here binding of Zac -repressor at the Operator 
sequence directiy destabilizes RNAP-promoter complex that in 
turn lead to the down regulation of transcription [1], [3]. The total 
time Xd^a required for the formation of a full-fledged w^-mer via 
3D diffusion-controlled collisions and subsequent binding with the 



cis-regulatory sites can be calculated as follows. 

(6) 

= naXs,\ + T^dA («<! + 1 - + 2) - yj 

In this equation riaXs^x is the time required for the searching and 
binding of the entire ?i„-mer at the corresponding CRMs on DNA 
(for a monomer it will be T, i) via a combination of ID and 3D 
diffusion, T^j is the time that is required for the formation of a 
dimer of the end-product through 3D diffusion under in vivo 
conditions, ^{na)=d\nr{na)/dn^ where T{na) is Gamma 
function and y^, = 0.5772157. is the Euler-Mascheroni constant. 
Here Td,i = 10'^ lenRDrNACz (s) is the minimum possible 3D 
diffusion controlled bimolecular collision time inside the cellular 
volume where R is average radius of the monomers of end- 
product, Cz (mol/lit) is the concentration of end-product 
inside the cellular volume, N^t is the Avogadro's number, 
Dt = kgT /6n<pK ( m^s ') is the 3D diffusion coefficient associated 
with the dynamics of monomers of end-product in aqueous 
medium where ks is the Boltzmann constant, T is the absolute 
temperature (K), (p is the viscosity of the medium. In the 
calculation of Trf_i, we have assumed that the reaction radius 
between two monomers is ~2R. Since the overall maximum 
radius of the m-mer will be ~inR, we find that Trf „,oc»i/(l 
and subsequently the total time that is required to form a Jio-mer 
in a sequential manner via 3D diffusion will be given by the sum 
Zd^i Ylm^l 'w/(l -|-m). We should note that this is the maximum 
possible search time since we have assumed a maximum possible 
radius for the n„-mer complex and also we have not considered the 
possibility of formation of the final Wa-mer through non-sequential 
and random pathways and the steric factor associated with the 
multimerization reaction. The total search-time that is required by 
the monomer of end-product to find its cognate site on DNA is 
defined as Ts,l =N{tl + t„s)/L where the overall ID sliding time is 
defined as Zl = iJ j 6xd [28] . In this calculation we have assumed 
that the end-product searches for its binding sites on DNA via a 
combination of ID and 3D diffusion controlled collision routes. 
Monomers of the end-product undergo at least N/L numbers of 
cycles of 3D diffusion mediated association that is followed by ID 
scanning and dissociation where N is the size of the genomic DNA 
(base-pairs, bps) and L (bps) is the average ID sliding-length 
between non-specific 3D association and dissociation. Here T/, is 
the time that is required by the monomers of end-product to scan 
L bps of the genomic DNA via ID diffusion along the DNA chain, 
Xd (bps^s ') is the ID diffusion coefficient associated with the 
sliding of monomers on DNA (this will be scaled down to Xd/n^ for 
w„-mer) and T„s is the time that is required for non-specific binding 
of end-products with the genomic DNA via 3D collisions under in 
vivo conditions. When all the monomers of end-product search 
for their binding sites on DNA in a parallel manner, then the total 
time [26] that is required {T:s,nJ for all these monomers to 
assemble at the sequentially located CM-acting elements can be 
derived from the theory of combinatorial binding of transcription 
factors [27-29] witii DNA as follows. 

'^s.na =N{zlK + t„s) /Lxz,_in'^ (7) 

From this equation we find that the ID scanning time increases 
with the number of monomers in a power law manner as 

where typical value of the exponent seems to be a2/5 [26]. From 
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Eqs (6) and (7) one can compute the following ratio. 



From the theory of site-specific DNA-protein interactions we 
find that (t^ i/Tsj) « 10^ for n„ = 1 [26-29] which suggests that 
0„^ > 10^ for aU values of Eqs (7) and (8) suggest the pathway 
I is more efficient than pathway II. This means that though the 
diffusion limited multimerization of the monomers of the end- 
product is not a reasonable assumption for large values of n^, 
direct assembly of the monomers of the end-products of TFs on 
the sequentially located cis-regulatory DNA binding sites via a 
combination of ID and 3D diffusion controlled routes can be still a 
reasonable assumption even for higher values of n„. In this context 
we can replace the combinatorial binding of numbers of TFs 
with CRMs of template DNA with a single step (ra„+l)* order 
reaction as given by the first equation in Eqs (4). One should note 
that unlike the prokaryotic systems, most of the eukaryotic 
promoters are activated through a combinatorial binding of 
several TFs at the corresponding CRMs [27]. This observation 
suggests that GO oscillator can be still a feasible model that can be 
used to generate limit-cycle oscillations in the cellular levels of 
negatively self-regulated TF proteins especially in eukaryotic 
systems. 

Steady-State analysis of Goodwin oscillator 

The system of Eqs (4) has a fixed point Pa = t]a which is a real 
solution of the following polynomial equation of the order (ria+l). 

tiuKHa + ina/i^a+ Ka))"")- naPa=0; Pa = i^+ffaKa/i^a + Ka)) (9) 



The steady state values of other concentration variables (Z^, Ma 
and Xa) can be calculated using the fixed point rja as follows. 

Zas = Va/i^a + K); ^as = " >1aPa) Mas = t]aPa 



Using the Jacobian matrix e\ aluated around the equilibrium 
point Pa = t]a, the linearized form of the system of Eqs (4) near 
this equilibrium point can be written as follows. 



/Xa\ 

Pa 
\Zal 



( -g 0 

-X/Wa -l/Wa 0 0 

0 1 -{\+CTa) CTaK 

8 0 X/Sa -c I 



Aaha\(Xa\ 
Pa 



(10) 



r = g + C+{\+aa) + \/Wa 

S = g{\+aa)+Cg + c{\+aa)+{g + C+{\+aa))/Wa 

-(^aKI^a-Aa^/Va 

t={{^+cy,)g+cg+c(\+a^)lwa+cg{\+aa)-{Aad{\+aa)+\lw^lva 

-<yaK(g+'^/Wa)/t.a 
U = Aa/Va£aWa + {cg{ 1 + <7a) -Aa3{ 1 -|- (Ta)/Va)/Wa - tTa^ag/WaSa 



From the Routh-Hurwitz criterion for a biquadratic polynomial 
[30] we find that the equilibrium point of the system Pa = rja will 
be stable only when all the following inecjuality conditions are true. 
In other words there may be limit-cycle oscillations around the 
steady state only when any of these inequality conditions is not 
true. 



r>0; {rs-t)>0; {rs-t)t-r^u>0;u>0 



(11) 



Here the first inequality condition r > 0 will be always true since 
tJa'xl/Pa ^'^'^ subsequently we find Aa > 0. The second inequality 
condition {rs — t)>0 can be reversed only when either i<0 or 
rs<t for s>0. When the third inequality in Eqs (11) is not true 
then the biquadratic polynomial can have a complex root such 
that >'Rc ± iyim with positive real-part (here >'Re > 0 yi^n > 0) that 
results in the generation of limit-cycle oscillations of the 
concentration of TF protein A and the period of such oscillations 
[29-32] wiU be given by Tp=2n/yi,a- This means that the period 
of GG oscillator can be modified by tuning the lifetime (y^,^) of the 
protein product of TF gene A (since T = y,,^,/) though the value of 
yim is a function of other Group I paramc-ters (-ui,,, Va, £a) which are 
in turn linearly depend on y^^. Since the HiU coefficient term n„ 
presents only in the coefficient terms 5, t and u, the third inequality 
in Eqs (11) can be reversed by increasing the value oin„ for any 
set of Group I, II and III parameters. This means that the 
parameter space that is required to generate oscillations can be 
expanded by increasing the value of n^. Inequality conditions 
given by Eqs (11) for a stable motion of the dynamical system of 
Eqs (4) can be directly derived from the following Routh table 
(RTqq) [30] corresponding to the biquadratic polynomial (PI). 



Rgg = 



Here we have defined various matrix elements as follows. 



1 s u Y'^ 

r t Q 

s-tjr M 0 

t-rul(s-tlr) 0 0 Y 



(12) 



Aa = riaHai^a + ««)(!- Pa^la)/ la 

The coefficients associated with the characteristic polynomial 
Y'* + rY^+sY^ + tY + u = 0 {PI) of the Jacobian matrix defined 
in Eqs (10) can be written as follows. 



When 7101 = 0 then the steady state solution will be either 
asymptotically stable or unstable depending on the values of real- 
parts. From Eqs (10-11) we find that the system will be 
inconsistent near the frxed point both at very large as well as small 
values of Group I type parameters as (va,Wa,Ea) ->0 or oo. This 
means that there exists a critical range of these parameters to 
generate limit-cycle oscillations in the ceUular levels of TF protein 
A. One should note that Group I parameters appear in the 
denominator of definitions of various coefficients of the charac- 
teristic polynomial (PI) which means that the period of oscillations 
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will increase proportionately with respect to an increase in these 
parameters since we find >'imOcl/(Va,H'a,Ea). Contrasting from 
Group I, there exist critical or threshold values of Group II type 
ordinary perturbation parami^tcrs {y_^i,a„,Ka) below which the 
oscillations occur. As in case of Group I type parameters, there 
exist a critical range of values of {Aa,fJ.„) Group III to generate 
limit-cycle oscillations. The critical value of Hill coefficient 
Ha = c^a for a given set of parameters can be iteratively calculated 
by numerically solving the third inequality in Eqs (11) at various 
values of When the dynamics of promoter state occupancy is 
much faster than the rate of change in the concentrations of other 
variables then one can set Va»0 and the system of Eqs (4) 
reduces to the following form. 

WadMjdx = nJ{n„ + Z„"") - Ma 

dPjdx = Ma-Pa~aa(Pa-XaZa) (13) 
&adZal dT = Pa-{la + Ka)Za 



lim,,^0 r,, = {n„{la + Kar/Pa)""''^ 

lim^^^oo >7<, = l/^a;lim„„^oo ria = l/Pa (17) 

Eqs (15-17) suggest that strong binding of the end-product 
(fig => 0) at the promoter of TF gene A is required along with the 
conditions such as (k^, ffa)—^t ^f^d (Aq, Wa, 6a) = 1 to decrease 
the required critical Hill coefficient towards the minimum possible 
value as cfta * 9. When there is an additional dimerization step 
as in Eqs (4) and (5) then the resulting characteristic polynomial 
of the Jacobian matrix will be of fifth order as 
Y^ + rY* + sY^ + tY^ + uY + m = 0 (PIII) and the Roudi criteri- 
on that is required to generate oscillations can be written as 
follows. 

K{tL-rK)-mL^<0;K = u~m/r;L = s-t/r 



Most of the earlier studies on GG oscillator consider Eqs (13) 
as the base model however with the conditions such that 
(o'fl,Ka) = 0. Using detailed numerical simulations we will show 
later that this assumption is reasonably invalid. The corresponding 
Jacobian matrix around the steady state Pa = rjg can be written as 
follows. 



dx\ 



Pa U 1 -(1+C7a) aaXa Pa 

,ZJ \ 0 l/Sa -{Xa + Ka)/eJ\ZaJ 



(14) 



Here we have defined A'^ = naPa{X — lia>la){^a + Ka) and the 
characteristic polynomial associated with this equation is 
Y^ + r'Y^+s'Y + t' = 0 (PII) where the coefficients are defined 
as follows. 

r'=l+<7a + {Xa + Ka)/ea+l/Wa 

s' = {Aa + Ka){Wa+l)/Waea + {l+<Ja)/Wa + Ka(Ja/Sa (15) 
t'=Aa/ WaBa + {Xa + Ka{l+<Ja))/Waea 



The physiological ranges of various parameters associated with 
Goodwin-Griffith oscillators are hsted in Table 1. 

One-to-one dual feedback oscillators 

One can extend these scaling ideas for one-to-one negative 
feedback oscillator or toggle switches (Figure 2B1) and repressi- 
lator models (Figure 2C1). In case of one-to-one negative 
feedback oscillator number of end-product molecules of TF 
gene A bind with the CM-regulatory elements associated with the 
promoter of TF gene B and subsequently down-regulates whereas 
«j number of end-product molecules of TF gene B down-regulate 
the promoter of TF protein A upon binding with the correspond- 
ing cis-regulatory elements (Figure 2B1). The set of difTerential 
rate equations associated with the two TFs one-to-one feedback 
system can be written in the dimensionless form as follows. 

ndXhldz = Z^'"i{\ -Xh)-HhXh 

WhdMiJdx = (\—Xh) — Ml, 

(18) 

PhdPh ldz = Mh-Ph-(yh[Ph- hZh) 

ShdZh/dz = Ph- {h + Kh)Zh - Xh {Zh"'' (l-Xg)- HgXg) 



The Routh-Hurwitz [30] condition required by the system of 
Eqs (13-14) to generate limit-cycle oscillations will be 
(/i' — t') < 0. Upon solving this inequality for the HiU coefficient 
Ha, the expression for the critical value of that is required to 
generate limit cycle oscillations can be obtained as follows. 

Cfla = {WaBar's' -(/.„ + ( 1 -|- <Ja)))/Pa{ 1 " Ma) + Ka) (16) 

Here one should note that the term )j„ in the right hand side of 
this equation is stiU a function of and fig and the following 
limiting conditions exist. 



In these equations the subscripts will be such that when h = a,b 
then q = h,a where (a, h) denote the TF genes A and B respectively. 
One should note that the HUl coefficients associated with the 
binding of the end-product of TF A at the promoter of TF B and 
end-product of TF gene B at the promoter of TF gene A are 
and til, respectively and in general na¥^ni,. Here we have defined 
various other dimensionless variables and parameters as follows. 

'^ = ypat\ Ph =Ph/Phs; Mf, = mh/mhs; Z^ =Zh/phs\ 
Xh = Xh/di,z;h = a,b 



<Th = Mf/yph ; Xh =p"hs ' dq,kgf / Xgf 
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Table 1. Simulation parameters used to integrate Eqs. (4-5) of Goodwin-Griffith osci 
using lad system (/io = 4) of £ coli. 


later model as constructed in reference [25] 




Parameter 


Definition 


Physiological values in £ coli 


Remarks 


Va 


ypa/pl'M 


-0.0001 


promoter state dynamics 


Wa 


Ifahma 


-0.5 


relative mRNA-protein lifetimes 


Ba 


ypa/Xaf 


-1 


end-product formation dynamics 


IJ-a 




-0.00002 


binding of end-product at promoter 


Xa 


Xar/Xaf 


-1 


end-product equilibrium dynamics 


Ga 


>-afh'fa 


-1 


connects protein and end-product formation 


Xa 




-1 


connects end-product and promoter state 

\Xy \ lai 1 11^^ 


Ka 


i-al ^af 


-0.01 


describes end-product decay dynamics 


"ya 


Pa,hyalK{ 


-0.1 






'ifalh-aPa, 


-0.5 




Xya 


Kya/PasXfya 


-1 






fi-r'da-M/Xfya 


-1 




daz 




1 


molecules 


Pa, 




-10000 


molecules 


f^as 


K,ah,„a 


-100 


molecules 


doi:l 0.1 371 /journal.pone.Ol 04328.t001 



hi = hr/ hf ; A'a = Khrf/pill ; Kh,-f = ky/khf, 



In these definitions for h = a, b one needs to set q = h, a. The 
steady state solutions {y\i,'^a) to the coupled set of Eqs (18) 
corresponding to this two TFs system witli respect to tlie scaled 
protein levels can be given as follows. 

P., = ^„ = e"; =nt= C'"' {p„ ( 1 - IUa)/nalhy/"' (19) 



The steady state values of other dynamical variables can be 
calculated using the steady state values of the protein products fj,, 
and rji, as follows. 

Zhs = ilJ {ki, + h);^hs = {i-tlhPh); Mi,s = ; h = a,h (20) 



Here y in Eqs (19) is the real root of 
y{nanb-l)-n,,ln<i+ln{^j;;li,-\\-l},ey))=0 where we 
have defined the function (I> as, 



Eqs (18) have three possible steady state solutions viz.()7j = )7„), 
{'lb<^la) 'ii'id >*?«)■ Under identical values of all the 
parameters such as (ji^ = jii„ n„ = ni, and so on) we find the 
unstable steady state solution of the two TFs system as rj/, = 
where 0< (j;^ = ly^,) < 1. This means that the limit-cycle oscillations 



around this unstable steady state can occur only when the values of 
all the control parameters and initial conditions are identical with 
respect to both the TF genes A and B. Using the eighth order 
characteristic polynomial of the Jacobian matrix associated with 
the linearized form of Eqs (18) (Methods section) near the steady 
state values (>76 = '?„), one can numerically derive the conditions 
for the occurrence of oscillations from the Routh-Hurwitz 
criterion. When the values of the control parameters are such 
that (fi„ ^ fij, or Ha # Wi and so on) or there is a transient 
perturbation in the values of these parameters or initial conditions, 
then the oscillating system wiU be unstable and driven to any one 
of the stable steady state solutions as either rji, < rj^ or rji, > r}^ 
through asymptotic spirals. For example when «/, » «„ or » 
then the stable steady-state solution will be (>7„«0, r]i,x 1). These 
results suggest that contrasting from GG oscillator model the 
identical two-TF feedback system cannot generate self-sustained 
oscillations in the presence of stochastic noise. One can also 
construct one-to-one feedback oscillator via coupKng two inde- 
pendent GG oscillators. Here these independent TF oscillators A 
and B can be coupled via either A-OR-B (Figure 2B2) or A- 
AND-B (Figure 2B3) type logics. One can consider various types 
of regulatory combinations associated with these network archi- 
tectures. The combinations in A-OR-B type coupling can be 
denoted as 'AsAc-BsBc' where 'As' and 'Bs' are the types of self- 
regulation of TF genes A and B respectively whereas 'Ac' and 'Be' 
are the types of their cross-regulation on each other. Each type of 
regulation can be either 'P' or 'N' where 'P' denotes positive type 
and 'N' denotes the negative type regulation. Using these notations 
one can denote the configuration given in Figure 2B2 as NN-NN 
type. Figure 2B4 as NN-PP type and Figure 2B5 as NP-NP 
type. The configurations given in Figures 2B4 and 2B5 are the 
well-studied robust dual-feedback oscillators [25], [33-35]. Sim- 
ilarly one can consider various possibilities in A-AND-B type 
architectures. Noting the symmetry of regulation we find three 
possible types as P-P, N-N and N-P out of which only N-N wiU be 
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a robust oscillator. The configuration given in Figure 2B3 is an 
N-N type. When the coupling is via A-OR-B type logic then both 
the promoters of TF genes A and B will be independently down- 
regulated upon binding of protein end-products of both the TF 
genes A and B at the respective cis-regulatory elements associated 
with each promoter and the first and last equations in Eqs (18) 
will be modified as follows. 



yhqdXhq/dx = Zg"'"l{\-Xh)—lli^^Xl^; Xh = X/7n = aA 

h,q=a,b 

ZhdZh/dx = Ph- {h + Kh)Zh- thq (Zh^ ( l-JTj) - ^hq^hq} 

-X,,{Z,"hh{l-X>,)-n,,X>,h) 



(21) 



In these equations for each value of subscript 'A' the subscript 
will take a, h and there are totally four equations associated with 
the (ncrall promoter state dynamics. Various modified parameters 
in Eqs (21) are defined as follows. 

A**? = Khqrf j Ph^ ; Kh^f = khqr /khqf;\hq= y^a j Phs ^hqf \ 

Ihq =Pq"' dq^khqf j Ihf 



The steady state solutions corresponding to the modified 
equations can be obtained by numerically solving the following 
set of coupled polynomial equations. 



l^hhl^hql {i^mP-hq + l^hhZcis"'"! + Hl,„Zh"''l' ) 

-hlh = ^'^h = a,b;q = b,c 



(22) 



Here Z/j, is defined as in Eqs (20). When fiij^ = fi and fii, = \ 
then we can calculate the steady state protein levels from the set of 
polynomial Eqs (22) as = (A/, + K/,)e^ where y is the real root of 
y(na + \) - ln(/i(l -e-)-) -e->'(«» + i)) =0. When the GG oscillators 
A and B are coupled through A-AND-B type logic then the dimer 
(y^) of both end-products (Za-Zj) wiU be the key regulating molecule 
that binds at the cis-acting elements associated with the promoters 
of both TF genes A and B however with different Hill coefficients 
(m/j) and subsequendy down-regulate them. The respective 
modified differential rate equations corresponding to the dimer- 
ization and binding of dimer at the promoters of TF genes A and B 
can be written as follows. 

ZddYa/dx = Z„Z,,-^j Y, - ^^^^^^ x,{ Y/h{\ - x^) - li.X,,) 
VhdXhldx= Y/i'(\-Xh)-HhXi,;h = a,b (23) 

EhdZh/dx = Ph- {h + Kh)Zh - XdhiZaZb - 4>d Irf) 



The modified and new parameters and variables in Eqs (23) 
are defined as follows. 



Yd =yd/Pas; Bd = Vpa/PbAdf, Xdh = ^dfPqs/ hf, 
'I'd = ^dr/ ^dfPbs'i Ih =Ph^ ^dfizkhf j ^hfPqs 



In the definition of //j for h = a,h one needs to substitute q = h,a. 
Here V (M-'s"') and Arf, (s"') are the forward and reverse rate 
constants associated with the diffusion limited dimerization 
reaction between the protein end-products of TF gene A and B. 
The corresponding steady state solutions to Eqs (24) can be 
written as follows. 

Yds = ZhsZqs j hq ; Zhs = Phsl {h-\-Kh); 

Xhs=YdJ'hl{nh+YdJ'h);Phs = nh (24) 



Here /j;, is the solution to the set of following polynomial 
equations. 

A*/,/ (a*/, + {h + Kh) (Aj +Kq))/ ^d)"'' ) ~ Phni, = 0; 

h = a,b;q = b,a (25) 



In this set of equations we need to set q = b for h = a and for 
A = ft we need to set q = a. The parameters associated with dual 
feedback oscillators are summarized in Table 2. 

Three gene repressilator type oscillators 

Similar to Eqs (18) one can write the set of differential rate 
equations associated with the three TFs repressilator model as 
follows (Figure 2C1). 



VhdXtf/dz = Zs"" {\—Xh)—[i^Xh 

WhdMiJdx = ( 1 — A/,) — M/, 

PhdPhjdx = Ml, - P;, - a/,{Pi,-Ai,Zh) 

ShdZh/dx = Ph- (aa -I- Kh)Zh - Xh {Zh"h (1 - JTj) - n^Xq 



(26) 



Here the subscripts wiU be such that 
{h = a,h,c;s = c,a,b;q = b,c,d) where (a, h, c) denotes respectively 
TF gene A, B and C in a cyclic {h, s, q) manner and the variables as 
well as various control parameters are defined as in case of Eqs 
(18) and generally «o Similar to Eqs (19) one can 

derive the steady state solutions with respect to the scaled protein 
levels for three TFs system as follows. 

Pas = ria = e^ 
Pbs = nb 



1 



-nc={l^aKt{^-kna)haPa) 



I lie 



(27) 



l/nc«6 



In this equations y is the real root of 
y{nanbnc + y)— ticnb in A + nc hi (j>— In \l/ = 0 where we have de- 
fined various other terms as, 
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Table 2. Various parameters associated with dual-feedback A-OR-B and A-AND-B type oscillators and their definitions. 



Parameter 


Definition 


Remarks 


V* 




Describes promoter state dynamics of TD gene 'h', for 
h = a, b, q = b, a 




>p/,/ i i'ih 


It^ldLIVt^llirtlMMpiL'Ldlllllt^Lllllt^^tJI IFUdlt; 11,11 — O, U 


eh 




end-product formation dynamics of TF gene h , h = a, 
b 


n, 




binding of end-product of TF gene A/B at promoter, 
h = a, b 


h 


hr/hf 


end-product equilibrium dynamics, h = a, b 




hi 1 7 pi, 


connects protein and end-product formation, h = a, b 




yd,/hf 


describes end-product decay dynamics 


h 


\+ai,Ki,l(h + Ki,) 




Sd 


ypa/Pbshf 


Describes dimerization reaction between TF proteins 
A and B 


Yj 


yj/p,,. 


Scaled concentration of Za-Zf, dimer. 


tik 


hlfPqs/hf 






Khj hfPhs 


Described dimerization equilibrium 


th 


pIu ^ dji^kjif I hfPiis 




\ihq 


f^lwf 1 Phi 


splits into four types of fij,^^ in A-OR-B type coupled 
oscillator 




ypiilPla '^li'lf 


I'/j splits into four types of v/nj in A-OR-B type coupled 
oscillator 


Xhq 


Pif.s' dtizkiuif 1 hf 




Khqrf 


klup/khqj 


Describes binding of end-product of TF gene 'q' at 
the promoter of TF gene 'h' 


h 


hr/ hf 


Describes end-product formation equilibrium of TF 
gene 'h' 


Note: subscript a 


denotes TF gene A and b denotes TF gene B. 





doi:l 0.1 371/journal.pone.01 04328.t002 

Using these steady state values of protein products );/,, one can 
write the steady state solutions to other dynamical variables can be 
written similar to Eqs (20) as follows. 

Zhs = ili,/{Kh + h); Xi„ = ( 1 - 17/, A,) ; Ml,, = r]f,lii,; h = a,b,c (28) 

Under identical values of all the control parameters such as 
{^^ = fif, = fic, na = ni, = nc and so on), the system reaches the steady 
state as {t]a = ^b — 1c) which is an unstable fixed point since even a 
small perturbation in the parameter values or initial conditions wiU 
drive the system towards a stable limit-cycle. As depicted in 
Figures 2C2 and 2C3 one can also construct the repressilator 
type model by cyclically coupling three independent GG 
oscillators A/B/C through -AND- or -OR- type logical gates as 
we have constructed in Figures 2B2-3. When the type of 
interaction is through -AND- type logic, then the z„-Zj dimer 
down-regulates TF gene B, Zj-z,. dimer down-regulates TF gene C 
and the z^-z^ dimer down-regulates TF gene A. The set of 



modified difTerential equations corresponding to the configuration 
that is given in Figure 2C3 can be written as follows. 

^hkdYhiJ dz = ZhZk- hk Yhk - <Pi,k ( Kk ( 1 - ^/t) - P-kXk) ; 

Yhk=yiik/Phs; h = a,b,c; k = h,c,a 

VhdXh/dz = Y"JI (l-Xh)- n^Xi, ; h = a,h,c; k = c,a,h 

(29) 

ei,dZi,/dT = Ph- (//, + k/,)Za 

~ Xliki^hZk — Hk Yhk) — Xhq {ZhZci — Xiui Yiiii) ; 

k = h,c,a; q = c,a,b 



The steady state solutions can be obtained by numerical 
methods from the following set of algebraic equations. 

Y,,ks=Zi.Zk/hk; Xi„ = y;I;J {m, + Yll) ■ 

/(„/(/</, + (z,,,Z/„/;.,,A,)"'') - ft,'//, =0 (30) 



In these equations for h = a, h, c one needs to set k=c, a, b and 
various new and modified parameters are defined as follows. 
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hk = hkr/ hkfPks', <Phk = ^hkfdzkP"h I hkfPks'-, 
^hk=ypa/Pkskhkf; 'Ihklq = hklqfPk/qs/ hf 

When the t^pe of interaction is through -OR- type logic as 
depicted in Figure 2C2 then the end products of both TF genes 
A and B can independently down-regulate TF gene B and the end 
products of TF genes B and C can independendy down-regulate 
TF gene C and so on. The set of modified differential equations 
associated with such system will be similar to that of Eqs (21-23) 
where the indices wiU be extended for three TF genes A/B/C. 
One can also consider a fuUy interconnected network of TF genes 
A/B/C. In these configurations the self-regulated promoters of 
each TF gene A/B/C will be negatively regulated by the end- 
products of the remaining two other TF genes. Here the mode of 
overall combinatorial interactions among these regulating end- 
products and the corresponding promoters can be either through - 
AND- or -OR- type logics as represented by the dashed lines in 
Figures 2C2-3. In case of -OR- type logical gate, the negatively 
self-regulated promoter of TF gene A will also have cis-regulatory 
binding sites for the end-products of both TF genes B and C and 
so on. One can write the modified set of differential equations 
associated with such fuUy interconnected configuration (Fig- 
ure 2C2) as follows. 

fhqdXi^/ d-c^Z^"' ( \- Xi,) - iii,,lXi„i: Xi, = X/,„ Ji,q,k = a,b,cEi,dZi,/dT 

_ fi'*-a*+Kj)z,-»,,(z;''''(i-z,)-ft,,,jf,,,)-x,,,,(z;'''*(i-x*)-^„jrM) 1(31) 

l-Z«(Z^"(l-J^/c)-M*A) j 

In the first one of Eqs (31), there will be three equations for 

each promoter and there are totally nine equations associated with 
the overall promoter state dynamics. In the second set of three 
equations as well as in the associatc-d parameters for each value the 
subscript 'h' from the set {a, b, c), the subscripts q and k will take 
the remaining values. This means that when h = a then q = b and 
k = r and so on. Various modified parameters in Eqs (31) are 
defined as follows. 

I^hq = Khqrf j P^^ I Khqrf =khqr/ hqf ; Vhq = j khqf 
Xhh =Plf ^dhzkhhf j hf] thq =p"q^ ^dqzhqf j hf, 

Xhk=plT'^dkzkhkf I hf 

The steady state solution to Eqs (31) needs to be obtained by 
numerically solving the following set of equations. 

I'-hhl^hql^hk I {j'-hhl^hql^hk + l^hhl^hq^Tl'^ + t^hht^hk'^ qs'' + Hhqt^hk^hs' ) 

-Phlh = 0; Z/i, = t]h/{h + Kh) 

In case of fully interconnected configuration through -AND- 
type logic that is depicted in Figure 2C3 (with dashed lines), the 
complex z^-Zb-Zc wiU be the key regulating molecule that binds with 
the promoters of all the three TF genes A/B/C and down-regulate 



them. Similar to Eqs (23) one can write the modified set of 
differential equations corresponding to repressUator configuration 
that is fuUy interconnected through -AND- type logic as follows. 

ZddY^ldx = Z,Z,Z,-<j>a Ya - Y,k=a,b,c ^^/Ml - ^k) - iikXk) 
VhdXhldx= Y/h(\-Xk)-n^Xh;h = a,b,c (32) 
&hdZh/d% = Ph- {h + Kh)Zi, - XjhiZaZbZc - Id Yd) 

Here we have defined Yj =ydlPm- The steady state solutions to 
this equation can be obtained by solving the following set of 
polynomial equations. 

i"/,/ {}^h + {ZhsZqsZks/ kd) "'')-Ph1h = 0; Z„s = nml {^m + Km) (33) 

Here various modified and new parameters are defined as 
follows. 

^d = Vpa/PcsPbshf, Xdh = hfPqsPksj hf\ 0d = hr j hfPcsPbs] 
Xh=P"as ' dhzkhf / hfPbsPcs 



In these equations similar to Eqs (31) for each value the 
subscript V from the set (a, b, c), the subscripts q and k will take 
the remaining values. This means that when h = a then q = b and 
k=c and so on for other values. 

Perturbation-responses of various gene oscillators 

Sample trajectories and phase portraits of GO oscillator for 
i'„ = (0,2 X 10""*) are shown in Figures 3A1 -3 and 4A1-3. 
Irrespective of the type of initial conditions and magnitude of 
the control parameters, the trajectories always start with an 
overshoot of protein production that is followed by asymptotic 
spirals towards a stable limit cycle. This seems to be an inherent 
property of negatively self-regulated loops [9]. Figures 3B1-4 
and 4B1-4 suggest that there exists an optimum range of Group I 
parameters Va = (4— 12) x lO"'^ and {wa,&a) e (0.2,1.8) at which 
the critical HiU coefficient (c«a) that is required to generate self- 
sustained oscillations is a minimum which is in turn strongly 
dependent on the promoter state occupancy parameter ji^. This 
optimum range is also dependent on the values of other Group II 
and III parameters. The optimum range of the conversion 
parameter seems to be h £ (0.6— 1.4). Results suggest that strong 
binding conditions ixjl X lO^'* (v^ 7^ 0) and < 10~ [^a = 0) are 
required to minimize the value of critical Hill coefficient with 
respect to changes in Group I type parameters. The minimum 
achievable values of critical Hill coefficients seems to be cWa = 6 
(Va#0) and c«a = 9 (Va = 0). When there is an additional 
dimerization step as described in Eqs (4-5) corresponding to 
the negative-feedback-only (NFO) model considered in reference 
[2.5], the minimum achievable critical HiU coefficient seems to be 
c«a = 3 . One should note that in the Lac I oscillatory system the 
effective Hill coefficient is Wa = 4 since four dimers of lac I end- 
products involved in the overall negative feedback. Numerical 
analysis of this NFO model system using the physiological range of 
parameters as given in Table 1 suggests that the period of 
oscillator can be well tuned by changing the promoter state affinity 
fi^ of the repressor without compromising the amplitude much as 
shown in Figure 3C1. 
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Figure 4B4 shows the strong influence of da on the critical cWa 
which means that the approximation (ffa = 0) as in case of most of 
the earlier studies on various genetic oscillators is not a valid one. 
At the critical Hill coefiicient, the period as well as amplitude of 
oscillations are strongly dependent on the Group I parameters as 
shown in Figures 5A1-3. These results also demonstrate how the 
oscillator responds to temporal perturbations is Group I param- 
eters. As we have shown in the theory section, the period of 
oscillations increases with increase in the Group I parameters 
whereas the amphtude seems to decrease as the value of Group I 
parameters increase. One should note that square of period of 
oscillation is inversely proportional to the total energy of an 
oscillator whereas the total energy is directly proportional to the 
square of amplitude. This means that the total energy of a GG 
osciUator can be fine-tuned by perturbing the Group I parameters. 
The Goodwin-Griffith oscillator sc(^ms to abruptly enter into the 
modified limit-cycle orbit upon introducing the perturbation and 
relax back much faster upon removal of perturbation in the 
parameter Vo = 0 rather than perturbations in other parameters 
{Wa,Ba). In the latter cases, as shown in Figures 5A2-3 the 
relaxation of oscillator to the original orbit upon removal of 
perturbation seems to be through slow asymptotic spirals. 
Figure 5B1 suggest that the period of oscillations increases 
monotonicaUy with respect to increase in the value of Group I 
parameters as we have predicted in the theory section. When 
Va #0 then there exists a range of Wq 6 (0.3,1) at which the period 
of limit cycle oscillations and the required c^a are almost 
independent of changes in w„. Figure 5B2 shows that when 
Va # 0 then the period of oscillations linearly increases as Ua 
increases whereas it linearly decreases with increase in Ca when 
v« = 0. 

The dual feedback motif (Figures 2B1-3) is an adaptable one 

that can act as a toggle switch as well as an oscillator depending on 
the type of configuration. As we have shown in the theory section, 
the configuration depicted in Figure 2B1 requires identical 
values of all the control parameters as well as initial conditions 
to generate coupled as well as synchronized oscillations. Particu- 
larly this configuration can efficiently act as a toggle switch since 
the fixed point i]„ = r\i, is an unstable one and even small 
perturbations in the parameters or initial conditions is enough 
for the system to exit from the synchronized limit cycle oscillations 
around this unstable fixed point and subsequently move towards 
any one of the two stable steady states. The configurations given in 
Figures 2B2-3 can act as coupled oscillators. Sample trajectories 
and phase portraits of one-to-one coupled oscillators correspond- 
ing to the configuration given in Figure 2B1 are shown in 
Figures 6A1-6. The minimum achievable value of the critical 
HiU coefficient that is required to generate self-sustained oscilla- 
tions around the unstable fixed point (f/^ = j/j) of dual feedback 
oscillator seems to be — 5 that is closer (^n;, = 6) to the critical 
value corresponding to GG oscillator. Variations of critical Hill 
coefficient with respect to changes in combination of different 
groups of control parameters are shown in Figures 6B1-8. 
Results suggest that the minimum value of critical HiU coefficient 
that is required to generate self-sustained oscillations can be 
achieved only when /i;,<2xl0~'' and (w/,,e/,) 6 (0.5,2). Fig- 
ure 6B1 also suggests that the inequality condition V/, >4 x 10^' 
is required for oscillations. From Figure 6B2 we find that c^h is 
also independent on the changes in the ordinary perturbation 
parameter Xa- However it is strongly dependent on and the 
condition (7/,<0.3 is required to achieve the minimum value of 
critical f-n/j. Results suggested that when there are identical 
perturbations in the given control parameter, then the one-to-one 
coupled oscillator (Figure 2B1) behaves similar to that of GG 



oscillator. That is to say the period of oscillations increases and 
amplitude decreases with an increase in Group I parameters. Here 
the identical perturbations are such that for the parameter Wh we 
have \vi, => W/, + i5,jv, where ii',, = W/, prior to perturbation and the 
magnitude of perturbation is such that dw^ = . When any of 
these two conditions fails, then the system wiU be driven towards 
the corresponding stable steady state. 

Upon receiving a transient pulse of perturbation or imbalance 
in the control parameters the dual feedback oscillator exits from 
the limit-cycle orbit with a time-delay (Tdei) and subsequentiy 
reaches one of the stable steady-states via asymptotic spirals. Here 
the target steady state is dependent on the type of disproportion in 
the parameter values among TF genes A and B. For example 
when the perturbation is from fi^ = fii, towards n^<Hi, then the 
target steady-state will be r\„ > since the binding of TF end- 
product B at the promoter of TF gene A is stronger than the 
binding of end-product A at the promoter of TF gene B. It seems 
that the value of this timodclay is dependent on the extent of 
disproportion (Ttt, where the subscript 'A' denotes the control 
parameter under consideration such as ii^j as well as duration of 
the perturbation (t,,,) in control parameters or initial conditions 
associated with the TF genes A and B. Here the percentage of 
disproportion or imbalance with respect to the parameter kh that is 
associated with TF gene B is defined as nk = \QQ\ka—kh\lki,. 
Figure 6C shows the variation of the time-delay with respect to 
changes in the extent of disproportion (7t^) in the control 
parameter fif^ and duration of perturbation x^. It seems that Xdel 
approaches zero independendy upon increase in both and Kk- 
Further simulation results suggest that the time-delay is 
independent on the time {Xpubr) at which the perturbation in the 
control parameter is introduced into the system. 

Contrasting from the configuration gi\ i;n in Figure 2B1, the 
limit-cycle orbits of the coupled oscillators depicted in Figur- 
es 2B2-3 are robust against transient imbalances in the control 
parameters. The minimum achievable value of the critical HiU 
coefficient seems to be cW/, = 4 for the osciUator with A-OR-B type 
logic (Figure 2B2) whereas (fih = 2 for the coupled oscfllators 
with A-AND-B type logic (Figure 2B3). Results suggest that the 
limit cycle orbit of coupk-d oscillators with A-AND-B and A-OR-B 
type logics arc stable one. Wh(;n there are temporal perturbations 
in Group I parameters associated with one of the Goodwin 
osciUators (TF gene A/B) then the other unperturbed oscillator 
responds to the changes in the behavior of the perturbed oscillator 
depending on the type of logical coupling between them. As shown 
in Figures 7A1-2, Bl-2 and Cl-2 in case of A-OR-B coupling 
an increase in the magnitude of Group I parameters associated 
with one of the osciUators A/B does not change the period of the 
entire system of oscillators (period-buffering) though there is a 
decrease in the amplitude of the osciUator that is perturbed in 
(w/,,e/j). The decrease in the amplitude might be partially owing to 
the period-buffering effect. In case of A-AND-B type logical 
coupling, increase in the magnitude of Group I parameters (ii'/,,e/,) 
increases the period of osciUations and decreases the amplitude of 
the entire system of osciUators that includes both TF genes A/B. 
Figures 8A1-4 suggest that an increase in the parameter Vh of 
one of the osciUators initially increases the amplitude of other 
oscillator to a maximum which then decreases later. Perturbations 
in Group I parameters (wv,,V/,) associated with one of the 
oscillators A/B also results in a phase-shift in cases of both A- 
AND-B and A-OR-B type logical couplings as shown in 
Figures 7-8A1-2 and Bl-2. Whereas perturbation in e/, affects 
only the amplitude and does not affect the phases of the coupled 
osciUators A and B as shown in Figures 7-8C1-2. Here one 
should note that in case of A-OR-B type coupling the parameter 
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Figure 3. Dynamical aspects of generalized Goodwin-Griffith oscillator. A1. Phase portraits of Goodwin-Grlfflth oscillator as described by 
Eqs (4). One needs to substitute Q = IVl (scaled concentration of mRNA) for red line, Q = X (promoter occupancy) for blue line and Q = Z (end- 
product) for pink line. Simulation settings are ju^ = 2 x 10"'', Va = 10~', {(Ja,Ka,la)=^ a"d we set {Ba,^a,Wa) = 1 which required a critical Hill coefficient 
of cflg = 6 to generate oscillations. Total simulation time is 1 00 (measured in terms of number of lifetimes of the protein product of TF gene A) and 
integration step is At=10~'. A2. Trajectories corresponding to the settings in A1. A3. Roots of the (biquadratic) characteristic polynomial (PI) 
associated with the Jacobian matrix for settings in A1. B1. Variation of critical Hill coefficient with the parameter set (/i„,Va). Minimum of this critical 
value seems to be achieved at ii^ ~ 10"*, and ~ 10~'. B2. Variation of critical Hill coefficient with the parameter set (/la^Ea)- With the optimized 
settings in B1, the system seems to be robust when e„ e (0.2,2). B3. Variation of critical Hill coefficient with the parameter set (/j^.w^). With the 
optimized settings in B1, the system seems to be robust when Wa e (0 2,2). 84. Variation of critical Hill coefficient with the parameter set (fia,^a)- 
Default values of other parameters in B1-4 are as in A1. CI. Variation of period and amplitude of the negative-feedback-only model considered in 
reference [25] with respect changes in the promoter affinity parameter fi^. Simulation settings are given in Table 1. Red solid line in the period and 
blue solid line is amplitude of the oscillator. 
doi:l 0.1 371 /journal.pone.01 04328.g003 



will be split into v/, => {vhh^'^hk) where we have the indices 
{h,k) =a,b as given in Eqs (21). Above results corresponding to 
A-OR-B type logical coupling with respect to changes in the 
parameter v/, are valid only when the temporal perturbations are 
the same for a given promoter of TF gene A/B. This means that 
for TF gene A (here we have subscript k = d) the extent of 
perturbation should be the same for both v^a and Vah while the set 
of parameters associated with the TF gene B (k = b) remains 
unperturbed. Here one should note that vm controls the dynamics 
associated with the binding of end-product of TF gene 'H' at its 
own promoter whereas Vhk controls the dynamics associated with 
the binding of end-product of TF gene 'K' at the promoter of TF 
gene 'H' as we have shown in Eqs (21). When there are 
perturbations in only one of these two split parameters (as 
Vu^(Vaa,Vu/,)) then the coupled system of oscillators seems to be 
dynamically unstable and also produces modulated beats as shown 
in Figures 7A3-4. The period of such beats increases as the 
imbalance in the set of split parameters v/jt increases as shown in 
Figures 7A5-6. These dynamical instabiUties as well as beats 
abrupdy disappear once the perturbations in Vhk are removed. 
Whereas the system of coupled oscillators relaxes back to the initial 
unperturbed limit-cycle orbit through asymptotic spirals upon 
removal of perturbations in case of Group I control parameters 
{wh,Eh). Results from Figures 7 and 8 suggest that coupled 
oscillators with -AND- type logic are more robust against 
promoter state perturbations than the -OR- type coupling. Period 
of a network of oscillators can be easily fme- tuned by manipulating 
merely one of the oscillators when the mode of coupling is via - 
AND- type. 

Sample trajectories and phase portraits of a repressHator 
configuration given by Figure 2C1 are shown in the Figures 
SlAl-5 in File SI. The minimum achievable value of the critical 
HiU coefiicient for the repressilator seems to be c^h = 2 similar to 
that of a one-to-one feedback oscillator with A-AND-B type logical 
coupling. As we have shown in the theory section, the steady state 
fixed point is more stable when the parameters and initial 
conditions are identical for all the TF genes A/B/C. When there is 
a transient perturbation in the control parameters then the system 
leaves the steady state and enters into a stable limit-cycle orbit 
through asymptotic spirals with a time delay Zjel as shown in 
Figure S1A4 in File SI. Here the magnitude of this time delay 
seems to be cfirecdy proportional to the extent of imbalance or 
disproportions in the parameter values as shown in Figure S1A4 
in File SI. Perturbation in the control parameter v/, associated 
with any one of the TF genes A/B/ C results in the decrease of 
amplitude of the perturbed as well as the one that regulates it. 
However perturbation in Vh does not affect the period of 
oscillations of the entire system of TF genes as shown in Figures 
SlBl-2 in File SI. This means that when is increased then the 
amplitudes of oscillations of TF genes A and C decrease whereas 
the amplitude of B is not affected. Perturbation in the control 



parameters (w/jje/j) associated with any one of the TF genes A/B/ 
C decreases the amplitude of oscillations of the TF gene that is 
regulated by the end-product of the perturbed TF gene and 
increases the amplitude as well as width of oscillations of the TF 
gene that is regulating the perturbc'd gx'ne. This means that when 
(WfljEa) increases then the amplitude of oscillations of TF genes A 
and B decreases whereas the width and ampHtude of TF gene C 
increases. Further results show that an increase in {wh,Sh) of any 
one of the oscillators increases the period of oscillations of the 
entire system of oscillators as shown in the supplemenl|^^ 
information Figures S1B3-4 in File SI. 

Sample trajectories and phase plane portraits associated with 
the configuration given in Figure 2C2 are shown in Figures 
S2A1-3 and SBl-4 in File SI. Contrasting from three TF genes 
repressilator model (Figure 2C1) the configurations given in 
Figures 2C2-3 do not require any asymmetry in the values of 
control parameters or initial conditions to trigger the stable 
oscillations. When the mode of coupling of TF genes A/B/ C of 
GG oscillators is through -OR- type logic then in the presence of 
identical values of all the sets of control parameters the TF genes 
A/B/C oscillate in a synchronized manner with respect to period 
and amplitude. When there is a perturbation in set of the control 
parameters Vhq (here v/, wiU be split into {vhh,Vhq) for each 
promoter) associated with any one of the TF genes A/B/ C then 
there are at least three difierent phases of responses. In the first 
phase, as shown in Figures S2B1-2 in File SI the system tries to 
resist the perturbation by keeping the synchronized limit-cycle 
orbit intact whereas in the second phase the system becomes 
unstable and chaotic whose magnitude depends on the extent of 
perturbation. Upon removal of perturbation, in the third phase the 
system enters into new asynchronous limit-cycle orbit with stable 
phase differences among the TF gene oscillators. When there is a 
perturbation in one of the split parameters (yf,h,Vhq), then as 
shown in Figure S2B1 in File SI the second phase will have 
several repeating elements of resistance and instability. Perturba- 
tion in the control parameters {WhM) seems to have similar efiects 
which are evident from Figures S2B3-4 in File SI. Contrasting 
from these results, the oscillator depicted in Figure 2C3 seems to 
be more robust against changes in the Group 1 control parameters 
and also they return back to the initial coherent type limit-cycle 
orbit upon removal of perturbations as shown in Figures S2C1-2 
in File SI. Sample trajectories of fuUy interconnected configu- 
rations given in Figures 2C2-3 (with dashed lines) are shown in 
Figures S3A1-3 and SBl-3 in File SI. These results suggests 
that the fully interconnected three TF genetic oscillator will be 
more stable against perturbations in the critical control parameters 
when the mode of coupling is through -AND- type logic than the - 
OR- type logic. 

Tuning capabihties of A-AND-B (Figure 2B3) and A-OR-B 
(Figure 2B2) type coupled dual feedback oscillators are demon- 
strated in Figures S4A-D and S4A-D in File SI. These results 
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Figure 4. Dynamical aspects of standard Goodwin-Griffith oscillator. A1. Phase portraits of the regular Goodwin oscillator as described by 
Eqs (13). Q = M for red line and Q = Z for pink line. Simulation settings are (o-„,iq,) = 0, /i„~10^'-and (c„,/i„,ii'„) = 1 which require a critical Hill 
coefficient of cno = 9- Total simulation time is 100 (measured in terms of number of lifetimes of the protein product of TF gene A) and integration 
scaled-time step is At = 10^^^. A2. Trajectories corresponding to the settings in A1 . A3. Roots of the (cubic) characteristic polynomial (Pll) associated 
with the Jacobian matrix for settings in A1. B1. Variation of critical Hill coefficient with the parameter set (/(„,ir„). With the optimized settings in B1, 
the system seems to be robust when ir,, e (0.3,2.5). B2. Variation of critical Hill coefficient with the parameter set (/(„,£„). With the optimized settings 
in B1, the system seems to be robust when £„ e (0.5,1.5). B3. Variation of critical Hill coefficient with the parameter set (/i„,/'.„). 84. Variation of 
critical Hill coefficient with the parameter set {K„,a„). Default values of other parameters in B1-4 are as in A1. 
doi:1 0.1 371 /journal.pone.01 04328.g004 



show that A-AND-B type coupled oscillators can be tuned by 
changing the promoter state binding parameter efficiently 
without conceding on the amplitude of oscillations. Increase in fii, 
monotonicaUy decreases the amplitude (and increases the period) 



of A-OR-B type coupled oscillators. Whereas A-AND-B type 
oscillator shows two distinct regions of responses with respect to 
changes in ^i/j namely a responsive region and nonresponsive 
region. For the settings in Figure S4B in File SI, the A-AND-B 
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Figure 5. Perturbation responses of Goodwin-Griffitli oscillator. A1 . Effects of perturbations in Va on the linnit-cycle orbit of GG oscillator. The 
default simulation settings are n^=2x\0~^ and Va = 5xl0~^, (<7a,Ka,Xa) = 0, {eaAa,Wa) = i which required a critical Hill coefficient of c"ci = 6- 
Perturbation introduced in the interval from time 30 to 1 00 by abruptly raising the value to = 10 x 10~'. Increase in Va increases the period and 
reduce the amplitude of oscillations. Total simulation time is 200 (measured in terms of number of lifetimes of the protein product of TF gene A) and 
integration step is At= 10~'. A2. Effects of perturbations in Wg. The default simulation settings are as in A1. Perturbation introduced in the interval 
from time 30 to 100 by abruptly raising the value to Wa = 2. Increase in Wa increases the period and reduce the amplitude. A3. Effects of perturbations 
in e^.The default simulation settings are as in Al. Perturbation introduced in the interval from time 30 to 100 by abruptly raising the value to Sa=2. 
Increase in Wa increases the period and reduce the amplitude. B1 . Effects of variation of on the period and critical Hill coefficient that is required to 
generate oscillations. With the current default settings, there exists a range of Wa e (0.5,2) at which the critical Hill coefficient is minimum. B2. Effects 
of variation of Ua on the period and critical Hill coefficient. Both period of oscillations and critical Hill coefficient are linearly dependent on o-^ and one 
cannot assume <7a=0 as in cases of several earlier studies. 
doi:l 0.1 371/journal.pone.01 04328.g005 



system responds to the perturbations in ^i, when lih<lO ^. When 
Hh>lO ^ then changes in fii, will not affect the period of 
oscillations much. On the other hand the amplitude (measured 
in terms of Ph/Phs) of A-AND-B oscillator is > 1 in a wide range of 
fih as well as Wh values. Further A-OR-B type coupled oscillators 
behaves similar to that of single GG oscillator with respect to 
changes in the perturbation of Group I parameter Wh- In case of 
both types of coupled oscillators the amplitude seems to be 
inversely proportional to the critical f-w,/, rather than the period of 
oscillations as it is apparent from Figures S4C and D in File S 1 
and there exists an optimum value of Group I parameter W/, at 
which the amplitude of oscillations is a maximum. Figures S5A- 
D in File SI shows that increase in Group I parameter w/, 
increases the period of oscillations and decreases the critical (Pf, in 
both -AND- and -OR- type coupled oscillators. However the 
period of -AND- type oscillator seems to be more robust against 
changes in than -OR- type coupled dual feedback oscillator. 
Contrasting from these there exist a cutoff value of Eii~2 (for the 
simulation settings in Figure S5 in File SI) beyond which the 
amplitude of oscillations is practically zero in both types of coupled 
oscillators that is apparent from Figures S5C and D in File SI. 
The discontinuities in plots depicted in Figures S4 and S5 in 
File SI mainly originate from the fact that (;ach data point 
belongs to different values of critical Hill coefficients. For example 
upon iterating the perturbation parameter, we first find out the 
critical Hill value for the perturbed values and then obtain its 
amplitude and period of the perturbed system with newly 
calculated hiU coefficient. This means that each data point in 
those plots comes from altogether different dynamical systems 
since the HiU coefficients are different. Therefore we may not 
expect continuous type plots either. 

Discussion 

Transcription factor gene oscillators play critical roles in driving 
cell-cycle to circadian rhythms. Here we have identified the critical 
control parameters associated with seff-sustained oscillations of 
such oscillators and classified them into Groups I, II and III 
depending on their functionality. Group I parameters control the 
intracellular dynamics of synthesis and degradation of various 
molecules associated with regulated TF gene (Figure 1). The 
parameter wi, of Group I describes the strength of coupling 
between the rate of degradation of mRNAs and the rate of 
degradation of corresponding TF proteins. We should note that 
transcription and translation of various TF genes of prokaryotes 
are simultaneously taking place well within the cytoplasm whereas 
in case of eukaryotes the transcription is taking place inside the 
nucleus and the synthesized pre-mRNA transcripts need to be 
spliced within the nucleus and then transported to cytoplasm after 
other post-transcriptional modifications through nuclear pores for 
translation. These differences in the cellular architecture demands 
higher Iffetimes (l/y^A) for eukaryotic mRNAs than the prokary- 



otic ones which results in the general observation that the values of 
Wh associated with various genes in prokaryotes are lower than 
eukaryotes genes [36], [37]. It seems that in case of yeast, the 
genome-wide \ alues of 117, varies from 0. 1 to 1 with a median of 
~0.3 [36]. These observations suggest that the naturally occurring 
or synthetic oscillatory motifs in the transcription networks should 
operate well within this range of Wh 6 (0.1,1) with different 
median values depending on the type of organism. For most of the 
bacterial genes we find that W/, ~ 0.1 [8], [9]. The parameter 
describes the dynamics of binding-unbinding of the end-product 
molecule with the promoter of TF gene A. Large values of Vq 
indicate slower changes in the residence state of promoters 
whereas small values indicate the faster dynamics of the promoter 
state towards the equilibrium. Earlier studies suggested [8], [9] a 
typical physiological value of as 10 ^ for a prokaryotic self- 
regulatory systems (w^ = 1) such as the one in E. coli. Whereas in 
case of eukaryotic systems such as yeast and human its value will 
be scaled up respectively to Vq ~ 10^^ and Va~ 10^^ owing to 
dilution in the local concentrations of various molecules upon 
increase in the nuclear volumes. One should note that yeast 
nucleus is ~10' times higher than E. coli cell whereas human cell 
nucleus is ~10^ times higher than E. coli cell. The parameter 
describes the dynamics associated with of conversion of the protein 
product to end-product towards the equilibrium. This conversion 
reaction can be either a first-order conformational transition or 
pseudo first-order chemical modification of the protein product via 
an additional catalyst. Here also acts as an indirect delay 
parameter that in turn relates the protein decay rate with the 
com ersion rate. One should note that the binding-parameter ji^^ is 
the central one that connects the entire Group I type singular 
perturbation parameters. By definition is inversely proportional 
to the affinity of the end-product towards the promoter sequence. 
Along with fi^ the ordinary parameters {('a,^a,Ka) decide the 
steady state value of TF protein. Further one should note that fi^ 
(or jU/, in general) is the only parameter that can be externally 
modified in a working model and all the others are fixed system 
parameters which can be modified only at the time of designing 
the oscillatory module. In case of Lac based oscillators ^la can be 
modified through changing the concentration of IPTG [25] which 
is a gratuitous inducer of lad gene. 

Understanding the design principles associated with the genetic 
oscillator is one of the central topics in synthetic and systems 
biology. An efficient synthetic oscillator should be robust against 
transient perturbations in the control parameters and fluctuations 
in the promoter state occupancies. The period of oscillators should 
be easily tunable without compromising the amplitude. Strieker 
et.al in reference [25] suggested that the robustness of the dual- 
feedback oscillators depicted in Figures 2B4-5 can be further 
increased by the explicit delay owing to the oligomerization steps 
associated with the end-products of TF genes A/B (here gene A 
can be lad and gene B can be araC). Further studies by Tsai et.al 
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Figure 6. Dynamics of one-to-one cross feedback oscillator. A1-2. Phase portraits of TF genes A and B which are coupled through one-to-one 
cross feedback loops. Simulation settings are {(T/,,k;,,7,,) = 0, {c/,,/t/,,»7,,p,,) = 1, = 10"^ and v;, = 5 x 10^^^ which required a critical Hill coefficient of 
cna = 6. Total simulation time is 100 (measured in terms of number of lifetimes of the protein product of TF gene A) and integration step is At = 10^^^. 
A3-4. Trajectories of TF genes A and B. Perturbation in was introduced at Tp,,/,., = 20 by abruptly raising ^i,^ to + 10^* for a period of t„ = 10^'. 
This corresponds to a disproportion of tt,, = 10^' (%). The system quits the limit-cycle orbit with a delay of t,/,/ ~50. Other default simulations settings 
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are as in A1-2. Q = X for red line, Q = IV\ for green line and Q = P for blue line. A5. Phase portraits of TF genes A and B. Q = P for red line, Q = IVl for 
blue line and Q = X for green line. A6. Roots of the eighth dimensional characteristic polynomial derived from the Jacobian matrix associated with 
Eqs (18) for the parameter settings given in A1-2. B1-8. Variation of critical Hill coefficient that is required to generate oscillations with respect to 
changes in various control parameters. Default settings of other parameters are as in Al -2. C. Variation of zjei with respect changes in percentage of 
disproportion tc^ and pulse width x„. 
doi:l 0.1 371/journal.pone.01 04328.g006 



in reference [35] suggested that unlike pure negative feedback 
systems, inclusion of a positive feedback loop within the negative 
feedback oscillators (Figure 2B4-5) can increase the robustness 
and tunable nature of the original oscillator without compromising 
on the overall amplitude. Here the positively auto regulated TF 
gene acts as a booster for the overall protein level of the negatively 
self-regulated oscillatory TF protein. In this context one can also 
consider a positive coupling of a positively auto-regulated TF 
protein K (booster) with N-N or NN-NN type dual feedback 
oscillators as depicted in Figure 2B6. Upon analyzing various 
classc-s of oscillators depicted in Figure 2 one can conclude that 
the effects of perturbations in the control parameters and 
promoter state fluctuations will be minimum when the mode of 
coupling among network of oscillators is through -AND- logic. 
From Figures S4 and S5 in File SI we fmd that the A-AND-B 
coupled dual feedback oscillator can be externally tuned by 
modifying the promoter state binding jii, without conceding on the 
amplitude of oscillations. The inherent chaotic response of the - 
OR- type coupled oscillators could be one of the reasons why the 
oscillations of most of the synthetic oscillators disappeared after 
certain period of time [18], [24] under in vivo conditions. 
Nevertheless the overall period of network of oscillators coupled 
through -OR- type logic are more resistant (period-buffering) to 
perturbations in the control parameters than the oscillators 
coupled through -AND- type logic. One should note that the 
temporal perturbations in the system parameters ultimately result 
in perturbations in the intra ceUular concentrations of protein end- 
products of coupled TF genes that in turn are looped I jack into the 
system. In this context the -OR- type coupled genes respond to 
perturbations in an additive way whereas those genes coupled 
through -AND- type logic respond in a multiplicative way. This 
could be the reason why the period of oscillations of -OR- type 
coupled oscillators are more resistant to perturbations than -AND- 
type coupled oscillators. This means that the period of network of 
oscillators coupled through -AND- type logic can be more easily 
tuned by perturbing any one of the coupled oscillators. Depending 
on the circuit functionality one can also chose a combination of 
both the types of coupling. So far we have assumed a copy number 
of TF genes as one (fi?^ =1). This may not be true in the bacterial 
based synthetic circuits constructed on plasmids since the plasmid 
copy number wiU be more than one in most of the times which can 
lead to further complexity. For example, a plasmid copy number 
of two for a NFO model that was constructed in reference [25] can 
mimic a combination of oscillators depicted in Figure 2B2 and 
2B3 with TF geru; A = B since it can mimic both -AND- as well 
as -OR- type coupled GG oscillators. This could be one of the 
possible reasons for observing stable and robust oscillations with 
such constructs [25] over several bacterial generations. 

In summary, in this paper we have considered various types of 
transcription factor oscillators by explicitly incorporating the 
promoter state dynamics and other chemical reaction balances in 
detail. Using our detailed model we have identified and classified 
various critical control parameters and numerically obtained their 
physiological and optimum ranges to generate self-sustained 
oscillations in the intracellular levels of mRNAs and transcription 
factor proteins. We further derived the basic design principles 
associated with robust and tunable gene oscillators. We have 



further demonstrated that by coupling two or more independent 
Goodwin-Griflith oscillators via -OR- or -AND- type logics one 
can construct genetic-osciUators which are fine-tunable and also 
robust against perturbations in the system parameters. When there 
is a perturbation in one of the -OR- type coupled oscillators, then 
the overall period of the system remains constant whereas in case 
of -AND- type of coupKng the overall period of the system moves 
towards the perturbed oscillator. Though there is a period- 
buffering, the oscillators coupled through -OR- type logic seems to 
be more sensitive to perturbations in the parameters associated 
with the promoter state dynamics than -AND- type. 

Materials and Methods 

We use Euler type numerical scheme [8], [9], [38] to integrate 
the set of differential rate equations corresponding to various types 
of oscillatory loops. For example in case of Eqs (4) which describe 
the Goorh\in-Griffith model, the numerical integration scheme 
can be written as follows. 

Xa,i+ 1 = Xa,i + [Z.f- ( 1 - Xa,d - llaXa,d^^lVa 
Ma,i+\ =^a,i + ((l -Xa^i) - M j W a 

Pa,i+ 1 = Pa.i + {Ma,i - Pa,i - <7a{Pa,i - ^aZa,i))Az (34) 
Za,i+1 — 

2a,i + {Pa,i - + - Zal^a,,"" ( 1 " X^,,) " /<„X«,,) ) At/e^ 

Initial and boundary conditions are (Xa,Ma,Pa,Za) =0 and 
{Xa,Ma,Pa,'Z.a) < 1 respectively. The scaled time step At should 
be chosen such that it captures the dynamics of all the variables 
including the dynamics of promoter state occupancies (first one in 
Eqs (34)). We divide the total scaled simulation time T into N 
equal intervals such that At = T/N. For simulation purpose we set 
At = 10~^ and the corresponding At = 0.03 s for a lifetime 1/ 
ypa~f>0 mins. We use Newton's method [38] to fmd the fixed point 
solutions to steady state equations. Here to obtain the solution to a 
nonlinear algebraic equation f{x) = 0 one uses the iterative 
scheme Xj+i =x,— /(x,)//'(x,). For example, the iterative 
numerical scheme to obtain the fixed point solution Pa = fla to 
Eqs (9) particularly for «a > 4 can be written as follows. 

'7a,; +1 

= fla,i+ [t^a- Pana,i[^ia + i^anZ))/{^^a + {^ta + lK';i^a)■, {^5) 

In this numerical scheme we set the initial guess value of the 
fixed point solution as )7^,-^q=0 and the tolerance limit for 
stopping iteration as Vla,i~''la,i+\\ — ^^~^- 
following computational workflow. For example, in case of one- 
to-one dual feedback oscillator one needs to construct the eighth 
dimensional Jacobian matrix associated with Eqs (18) as follows. 
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Figure 7. Dynamics of dual feedbacl< Goodwin-Griffith oscillators coupled through OR gate. A1-2. Trajectories of protein-products of TF 
genes A and B which are two independent GG oscillators coupled through A-OR-B type logic as given in Figure 2B2. Simulation settings are 

(f<;,A/i,Z;,,) =0, {Ri,,h,Wh,Ph) = 1' /ift = 10~'and Vh = 10~'which required a critical Hill coefficient of cna = 8. Total simulation time is 200 (number of 
lifetimes of the protein product of TF gene A) and integration step is z(t=10"'. For each promoter A/B the parameter f/, will be split into 
v/i=>(v/,;,,F/,^) where v/,/, corresponds to self-regulation and i'/„, corresponds to cross regulation. Under identical values of all the parameters the 
system generates synchronized oscillations with a period of t,, ~4.5. Upon introduction of perturbation in v„ from scaled time 0 to 100 (where 
Va = 15 X 10^^), the amplitude of TF gene A is reduced with a phase shift and the period of entire system that includes both TF genes A and B remains 
the same. A2 is a magnification of certain range of Al . A3-4. Effect of perturbation in only one of the split parameters (v/,/,,v/,,) associated with TF 
gens A/B. Here Vaa is perturbed to Vaa = 15 x 10^'. The system seems to be unstable and generates beats. A4 is a magnification of certain range of A3. 
A5-6. Here v„a is perturbed to Vaa = 30 x 10"\ Period of beats seems to increase as the disproportion among the split parameters increases. A6 is a 
magnification of certain range of A5. B1-2. Effect of perturbation in the parameter h'„ which is raised to = 1.5 from the default value in the interval 
from 0 to 100. Increase in Wa reduces the amplitude of both the TF genes A and B with a phase shift and the period of oscillations of the entire system 
remains the same as in A1-2. B2 is a magnification of certain range of B1. CI -2. Effect of perturbations in the parameter which is raised to Sa = 1-2 
from the default value in the time interval from 0 to 100. Increase in reduces the amplitude of both the TF genes A and B without a phase shift and 
the period of oscillations of the entire system remains the same as in A1-2. C2 is a magnification of certain range of CI. 
doi:l 0.1 371/journal.pone.01 04328.g007 
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Here subscripts [a, h) denotes the TF genes A and B respectively. 
Various matrix elements are defined as follows. 

CH = KH + h + Xh^k 



find that khf ~ 10^' molecules ' (y^) that wiU be scaled down to 
the interaction of W/, number of TF molecules with sequentially 
located combinatorial binding sites as k/,/ {khf /nh) 
molecules """S"' [26], [27]. The time associated with the 
unbinding reaction wiU be closer to that of the protein lifetime. 
Here we have assumed an in vivo 3D difiiision controlled collision 
rate ~10*M" s"'. In case of nucleus of yeast cell khf ~ 10 
molecules"' Yp~' and in case of nucleus of human cell we find 
khf ~ W^'^ molecules"' Yp~'- Since li„~10 * for a typical 
bacterial promoter [8], [9] we assume {ypi,/p"l;ki,f ) ~ 10"'* for 
an arbitrary Hill coefficient n„. Using these values one can 
estimate the physiological ranges of various groups of control 
parameters as given in Table 1. With these settings for a bacterial 
lad (Wo = 4) system (Eqs. (4-5)) the values of parameters will be 
Va<10^'* and /i^ ~ 10^'. For simplification purpose we can 
assume identical values for similar group of parameters Uh = Hhq, 
fill = and so on for other family of parameters such as Xh ^'^'^ 
Sh- 



in these definitions for h — a,b one needs to substitute k = h, a. 
To obtain the critical value of the Hill coefficient (cWa) one need to 
first solve the steady state equations. Upon substituting the steady 
state protein values into to the corresponding Jacobian matrix one 
can construct the characteristic polynomial and the Routh table. 
The corresponding inequality conditions for oscillations will be 
derived from this table. This procedure need to be carried out at 
various values of from = 1 in an iterative way. We set the 
default initial value of promoter state variable as X/; = 5 x 10"^ for 
any one of the TF genes A/B/C to trigger the limit-cycle 
oscillations in case of repressilator so that the effect of other Group 
I control parameters on the initial delay in oscillations can be 
studied. One should note that this initial delay is still a function of 
the disproportion in the initial conditions. We measure the 
concentrations of various molecules in terms of number of 
molecules inside the cell. Considering a typical E. coli bacterial 
cell (volume ~10~'" m'') [39] we set fi?/,^ = 1, m/j,~10^ molecules 
andj!>fo~10* molecules [8], [9] where h = a, b and c respectively 
denotes TF genes A, B and C. Concentration of a single TF 
molecule inside a bacterial cell will be ~2 nM [40-42]. We 
assume a lifetime of TF protein as ~2 min and mRNA lifetime as 
~0.2 min for E. coli that gives a decay rate of y^^ ~ 10^^.?^' so 
that = 0. 1 and we measure the time in terms of number of 
protein lifetimes. Since the dynamics of binding-unbinding of a 
single transcription factor protein {n„ = 1) with its CM-regulatory 
site on DNA is a typical diffusion-controlled site-specific DNA- 
protein interaction, under in vivo conditions of a bacterial cell we 



Supporting Information 

File SI This file contains Figure SI -Figure S5. Figure SI. 
Three gene repressilator model. A 1-4. Phase portraits and 
trajectories of TF genes A, B and C of a repressilator. Simulation 
settings are {oh,Kh,Xh) =Q, {shAi,Wh,Ph) = i, /^/, = 10"'* and 
V/, = 10^** which required a critical Hill coefficient of f;n„ = 2. 
Total simulation time is 200 (number of lifetimes of the protein 
product of TF gene A) and integration step is At=10~^. To 
trigger the oscillations, we have introduced the asymmetry in the 
initial condition for the promoter state occupancy of TF gene A as 
Xa = 5xlO^^. OscUlatioiis starts with a time delay Zjei whose 
value depends of the magnitude of this disproportion in the 
parameter values. A5. Roots of the twelfth degree characteristic 
polynomial associated with the Jacobian matrix of Eqs (26) for 
settings given in Al. Bl-2. Effects of perturbation in Va that is 
raised to i„ = 10^' (^a = 10^^ in B2) in the time interval from 0 to 
100. Increase in increases the period of oscillations of the entire 
system from Tp~23 to 24.5 and reduces the amplitudes of TF 
genes A and C. The amplitudes of TF genes A/B/C are such that 
A<C<B. B3-4. Effects of perturbation in {Wa,Ea) which are raised 
to It'a = 3 in th(" time interval from 0 to 100. Increase in Wa 
increases the period of oscillation of the entire system from Zp ~ 23 
to 30 and reduces the amplitudes of TF genes A and B and 
increases the amplitude of C and the amphtudes of TF genes are 
such that B<A<C (B3). Increase in Ea increases the period of 
oscillation of the entire system as in B 3 where the amplitudes of TF 
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Figure 8. Dynamics of dual feedbacit Goodwin-Griffith oscillators coupled through AND gate. A1-4. Trajectories of protein-products of 
TF genes A, B which are two independent GG oscillators coupled through A-AND-B type logic as given in Figure 2B3. Simulation settings are 
{oi,,Ki,,xi„j) =0, {Ei,,Xi,,wi,,ij>j,pi,) = 1, /I/, = lO^^^and v/, = 10"^ which required c''o = 2 (here we have set nj, = 4 for clarity of results). Total simulation 
time is 200 (number of lifetimes of the protein product of TF gene A) and integration step is Ar= 10~^. Under identical values of all the parameters 
the system generates synchronized oscillations with a period of t,, ~9. Upon introduction of perturbation in v„ from scaled time 0 to 100 (where 
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Va = 15x 10 ' In A1, Va = 30x 10 ' in A3), the amplitude of TF gene A is reduced with a phase shift and the period of entire system (bothTF genes A 
and B) increases to T;,~ 10 in A1. A2 and A4 are magnifications of certain range of A1 and A3. B1-2. Effect of perturbation in the parameter which 
is raised to =2 from the default value in the time interval from 0 to 1 00. Increase in Wa increases the period of oscillations of the entire system to 
Tp~9.5 as in A1-2 and reduces the amplitude of TF genes A with a phase shift. B2 is a magnification of certain range of B1. CI -2. Effect of 
perturbations in the parameter Sa which is raised to ea = 10 from the default value in the time interval from 0 to 100. Increase in increases the 
period of oscillations of the entire system to ~9.5 as in A1-2 and reduces the amplitude of both TF genes A and B without a phase shift. C2 is a 
magnification of certain range of CI . 
doi:1 0.1 371/journal.pone.01 04328.g008 



genes A/B/C are such that B<A<C (B4). Figure S2. 
Dynamics of three independent Goodwin-Griffith oscil- 
lators cyclically coupled. Al-3. Phase portraits of TF genes A/ 
B/C which are three independent GG oscillators cyclically 
coupled through -OR- t^-pe logic as given in Figure 2C2 (without 
dashed lines). Simulation settings are ((7/,,K/,,j;;,) =0, 
{eh,Xh,Wh,Ph) = i, = 3^11^ Vf, = lO^* which required a 

critical Hill coefficient of = 5 (we have set this to 6 for clarity of 
results). Total simulation time is 500 (number of lifetimes of the 
protein product of TF gene A) and integration step is At=10^^. 
In this configuration the end-products of TF gene A and C wUl 
regulate TF genes A through A-OR-C type logic whereas the 
promoter of TF gene B will be regulated by the end-products of 
TF genes A and B through A-OR-B type logic and so on. Under 
identical values of all the parameters the system generates 
synchronized oscillations with a period of T^~7. Bl-2. Effects of 
perturbation in v/,. Upon introduction of perturbation in Vaa from 
scaled time 100 to 400 (where Va = 6 x 10"'') there three phases of 
responses. In the first phase, the system tries to resist the 
perturbation whereas the second phase consists of repeating 
elements of resistance and chaos. Upon removal of perturbation 
the system enters into new limit-cycle orbit in the third phase. In 
B2 both Vaa and Vac are perturbed as in Bl. B3-4. Effects of 
perturbation in (vi^a>£a) which are raised to {ea,Wa)=2 in the time 
interval from 100 to 400. System responds to the perturbation as 
in Bl-2. Cl-2. Trajectories of TF genes A/B/C which are three 
independent GG oscillators cychcaUy coupled through -AND- type 
logic as given in Figure 2C3 (without dashed lines) and their 
responses to perturbations in Group I control parameters. 
Simulation settings are (ffA.KA,/^) =0. = 1, 

/i;, = 10^^ and vi, = \0^* which required cW,, = 2. Total simulation 
time is 500 and integration step is At = 10^^. Identical values of aU 
the parameters of the system generate synchronized oscillations. 
Introduction of perturbation in Va from scaled time 100 to 400 
(where Va = 3xlO-'', CI) affects only TF gene A whereas the 
orbit of other TF genes B/C seems to be stable and upon removal 
of the perturbation the system returns back to initial limit-cycle 
orbit. In C2 the parameter Wa is perturbed to Wa = 2 as in CI. 
Figure S3. Dynamics of three Goodwin-GrifBth oscilla- 
tors which are fully interconnected. Al-3. Phase portraits of 
TF genes A/B/C which are three independent GG oscillators 
which are fuUy interconnected with -OR- type logic as given in 
Figure 2C2 (with dashed lines). Simulation settings are 
(cta,Ka,j:a)=0, (Ei,M,Wh,Ph) = '^, /^A = 10"^ and v/, = 10"* which 
required a critical HiU coefficient of = 6. Total simulation time 
is 500 (number of lifetimes of the protein product of TF gene A) 
and integration step is At= 10~^. In this configuration all the end- 
products of TF gene A/B/C will regulate all the three TFs 
through A-OR-B-OR-C type logic. Identical values of all the 
parameters of the system generate synchronized oscillations. 
Perturbation in Va from scaled time 100 to 300 (where 
Va = 5 x 10~^, Al) seems to make the system unstable. In A2 Wa 



is perturbed to w a — 3 as in Al and m A3 Sa is perturbed to Sa — ^ 
as in Al. Bl-3. Phase portraits of TF genes A/B/C which are 
three independent GG oscillators fully interconnected with -AND- 
type logic as given in Figure 2C3 (with dashed lines). Simulation 
settings are (ffA,K/i,Xd/i) =0, {sh,^h,Wh,Ph) = l, /i/; = 10"^ and 
v/, = 10~^ which required a critical HiU coefficient of cWa = 2. 
Total simulation time is 500 (number of lifetimes of the protein 
product of TF gene A) and integration step is At=10^^. In this 
configuration all the end-products of TF gene A/B/C will be 
regulated by their complex. Identical values of all the parameters 
of the system generate synchronized oscillations. Introduction of 
perturbation in v„ from scaled time lOO to 300 (where v^, = lO^'' in 
Bl) seems to make the system unstable. In B2 M'a is perturbed to 
= 3 as in B 1 . In B3 Sa is perturbed to = 3 as in B 1 . Figure 
S4. Tuning dynamics of A-OR-B and A-AND-B types of 
coupled oscillators with respect to perturbations in /i 
and u>. A, C. Tuning capability of GG oscillators coupled 
through A-OR-B type logic as given in Figure 2B2. Default 
Simulation settings are (^ai„Ki,,Xhq) =0, {ei,,Ai,,Wij,piJ = 1 and 
v/, = 10^'. B, D. Tuning capability of GG oscillators coupled 
through A-AND-B t^pe logic as given in Figure 2B3. Default 
Simulation settings are ((T/,,K/,,X/,^) =0, (g/,,A/,,H'/,,^rf,P/,) = 1 and 
Va = 10^^. Plots A and B show the variation of period, critical cW/, 
and amplitude with respect to changes in jX/, (iterated from 
5x10 'to 10 ^ with Wi,= 1) whereas plots C and D show the 
variation of these quantities with respect to changes in W/, (iterated 
from 0.1 to 10 with ^;, = 10^^). Here period of oscillator is 
measured in the number of lifetimes of TF protein A (l/y^) and 
amplitude is measured in terms of number of Ph/Phs- Figure S5. 
Tuning dynamics of A-OR-B and A-AND-B types of 
coupled oscillators with respect to perturbations in v 
and e. A, C. Tuning capability of GG oscillators coupled through 
A-OR-B type logic as given in Figure 2B2. Default Simulation 
settings are ((7a,ka,Xa?) =0. (sa.Aa.wa.Pa) = 1 and ^</, = 10"^ B, 
D. Tuning capability of GG oscillators coupled through A-AND-B 
type logic as given in Figure 2B3. Default Simulation settings are 
{ch,Kh,Xhg) =0, {ehM,Wh,<l>d,Ph) = 1 and = 10~^ Plots A and B 
show the variation of period, critical f;W/, and amplitude with 
respect to changes in Vh (iterated from 5x10 ' to 10 with S/j = 1) 
whereas plots C and D show the variation of these quantities with 
respect to changes in Sf, (iterated from 0.7 to 8 with va = 10^^). 
Here period of oscillator is measured in the number of lifetimes of 
TF protein A (1/ya) and amplitude is measured in terms of 
number oi P^/Phs- 
(ZIP) 
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